Volcanic arc rigidity variations illuminated by coseismic deformation of the 2011 Tohoku-oki M9

Rock strength has long been linked to lithospheric deformation and seismicity. However, independent constraints on the related elastic heterogeneity are missing, yet could provide key information for solid Earth dynamics. Using coseismic Global Navigation Satellite Systems (GNSS) data for the 2011 M9 Tohoku-oki earthquake in Japan, we apply an inverse method to infer elastic structure and fault slip simultaneously. We find compliant material beneath the volcanic arc and in the mantle wedge within the partial melt generation zone inferred to lie above ~100 km slab depth. We also identify low-rigidity material closer to the trench matching seismicity patterns, likely associated with accretionary wedge structure. Along with traditional seismic and electromagnetic methods, our approach opens up avenues for multiphysics inversions. Those have the potential to advance earthquake and volcano science, and in particular once expanded to InSAR type constraints, may lead to a better understanding of transient lithospheric deformation across scales.


INTRODUCTION
Large-magnitude subduction zone earthquakes cause widespread deformation, including coseismic subsidence beneath volcanoes located hundreds of kilometers away from the rupture area (1,2).Associated local deformation may indicate rheological heterogeneity associated with the magmatic plumbing system (3,4) where the presence of high geothermal gradients (3,5) and crustal intrusions (6,7), for example, may mechanically weaken the rocks.
While geodetic observations have previously been used for inferences about elastic structure on global scales (8), here, we present the first direct inversion of coseismic crustal deformation for both upper mantle, margin-scale elastic properties and fault slip.We focus on northeastern Japan and invert coseismic surface displacements due to the 2011 Tohoku-oki M9 earthquake to investigate rigidity variations along the margin and specifically within the overriding plate's volcanic arc.Low-seismic velocity beneath the volcanic arc in northeastern Japan have long been imaged by seismic tomography (9,10), including for individual volcanoes such as Naruko (11,12), and high v P /v s ratios at depth suggest the presence of fluids and/or partial melt (10,13).The sources of these fluids may be attributed to dehydration metamorphic reactions of minerals from the subducting slab at ~100 km depth, above which many of the volcanic arcs worldwide are typically located (14,15).
Beneath the volcanoes and inland northeastern Japan, the upperplate structure exhibits shallow seismicity within the brittle deformation regime.In addition, deeper low-frequency earthquakes occur around active volcanoes, often associated with magmatic activity (9).In contrast, seismic activity is notably elevated and more widespread offshore, in particular along the primary thrust region beneath the landward slope of the Japan Trench.Several studies involving seismic, gravity, and residual topography analyses have linked the segmentation of megathrust events to material heterogeneities in the overriding plate (16,17), suggesting a substantial influence of structural heterogeneities on earthquake nucleation and rupture style (18).For instance, Bassett et al. (16) proposed that variations in the forearc lithology of the upper plate played a key role in the along-strike variations in the 2011 Tohoku-oki M9 earthquake slip distribution.It is thus important to independently quantify the extent and distribution of the Earth's rigidity structure within a fault and volcano dynamics context.

Joint inversion of geodetic data
We seek to infer variations in elastic parameters directly from coseismic geodetic deformation observations (19) associated with the 2011 Tohoku-oki M9.We use three-dimensional displacement information from 1283 permanent Global Positioning System (GPS) sensors, eight seafloor acoustic GPS sensors (GPS/A), and six pressure gauges (APG); the latter measure verticals only (20)(21)(22).The M9 resulted in major deformation (20)(21)(22), with ∼5m displacements recorded horizontally in the eastern part of the Tohoku region on land, and over 30 m on the seafloor (Fig. 1A) with data corrected to remove aftershock effects (23,24).Vertical deformation showed ∼5 m of uplift near the Japan Trench and ∼1 m of subsidence near the Pacific coast of Honshu.Observational errors for the inland and offshore stations are inferred to range from a few millimeters to several centimeters (table S1).
We perform an inversion of the geodetic data; first for fault slip alone, and then for slip and three-dimensional (3D) rigidity variations (Materials and Methods).Our model domain covers a 3700 × 4600 km rectangular region around Honshu with a depth of 700 km.We refine the mesh around the epicenter and the central part of eastern Honshu (fig.S1), where we anticipate the highest information recovery from the geodetic data (Materials and Methods).Figure 1B shows the inferred fault slip and residuals between model predictions and horizontal and vertical coseismic deformation of our best homogeneous structure model.The misfit is at the ≈2 and ≈26% level for the horizontals and verticals on land, respectively (Fig. 1B).The horizontal and, more so, the vertical residuals are clearly not random but show spatially coherent patterns.Such residuals are also apparent in earlier homogeneous medium inversions (25).Their origin in terms of a contrast between strong slab and weak overriding plate and/or mantle has previously been explored by forward computations exploring the role of elastic heterogeneity (24)(25)(26)(27).
Figure 1C shows that the residuals can be markedly reduced if we allow for 3D rigidity variations during a joint inversion for fault slip and rigidity, with the corresponding shear modulus variations shown in Figs. 2 and 3.In this case, root mean square (RMS) misfit values between model and geodetic data on land are roughly halved, to ≈1 and ≈16% for the horizontals and verticals, respectively (Fig. 1C and table S1).By jointly inverting for rigidity variations and coseismic slip, the remaining misfit patterns now lack any obvious, overall coherent structure.Only a small region in the eastern part of the Onikobe volcanic area shows overestimated vertical displacements.
The residual reduction is possible because deep elastic modulus variations have subtle but robust diagnostic surface deformation signatures for a given slip geometry (19,27), and the data appear to resolve an overriding plate that is weaker than the slab and additional anomalies in the mantle wedge that are associated with the volcanic arc (Fig. 2).Comparing the slip distribution inferred assuming a homogeneous structure with the joint inversion (contours in Figs. 1, B and C, and colored figs.S2, A and B), we find similar slip distributions with most slip concentrated near the trench, consistent with previous inversions with prescribed 3D shear modulus variations (24,26).The preferred joint inversion estimates a slightly lower maximum slip of 41 m compared to 46 m in the homogeneous case and has relatively smoother slip distribution, highlighting the trade-offs with fault-proximal structure.This reduction in slip magnitude broadly aligns with previous studies considering heterogeneous material structure in forward tests (24,26,27).Higher slip near the trench might be likely, e.g., from bathymetric surveys (28), within uncertainties (29,30), for example, and additional constraints could be included in future slip inversions.In our tests, we also find a moderate trade-off of maximum slip with slab geometry (fault dip; figs.S2C and S8), whereas the 3D rigidity structure is quite stable with respect to fault dip (figs.S6 and S7).

Rigidity structure and weaker volcanic arc
Concentrations of more compliant material underneath the Japanese landmass correlate well with the location of volcanic centers between 37° and 42°N (Fig. 2) where GPS spatial coverage and resolution are high (fig.S3).On the basis of a number of synthetic tests, we found these features to be robust outcome of the inversion (e.g., figs.S4 and S5).At shallow depths (Fig. 2A), rigidity reductions of ∼15 to 25% are observed around Mt. Akitakoma (40°N), Mt.Kurikoma (39°N), Mt.Zao (38°N), and reductions between 25 and 35% are seen beneath Mt.Azuma (∼37.5°N) and Mt.Nasu (37°N).These volcanoes experienced notable subsidence during the 2011 Tohoku-oki event (1).The presence of high geothermal gradients (3,5), and presence of deep-seated hot plutonic bodies (6) in this area aligns with the rigidity reduction due to increased temperatures weakening the rocks.These weaker material anomalies are prominent within the uppermost 30 km, which corresponds to the average crustal thickness in northeastern Japan (31) and to the region of reported shallow seismicity (black dots in Fig. 3), some of which may be associated with the volcanic plumbing systems.Seismic tomography of the Japan arc (32,33) suggests a correlation between the more compliant anomalies we image and low seismic velocity anomalies at shallow depths (Fig. 4).
An interesting observation in profile C-D of Fig. 3 is a lower rigidity region with a ~30% shear modulus reduction between 30 and 50 km depth beneath Naruko volcano, surrounded by a slightly stronger region (reduction below 20%).This deeper anomaly is well correlated with low-velocity anomalies and high v P /v s ratios observed beneath this volcano at these depths (11,12,34), possibly indicating the presence of a deeper magma reservoir.This deeper magma storage beneath Naruko, as well as shallower reservoirs beneath other volcanoes, may be directly connected to the source region of melt and fluids associated with metamorphic reactions within the subducting slab (35).
Notably, the volcanic low-rigidity anomalies in Figs.2B and 3 align with the subducting slab interface reaching ~100 km depth, and the equivalent distance from the trench is typically associated with the main volcanic arc (14,15).Shear modulus reductions between 8 and 12% are apparent at 75 km depth (Fig. 2) along the volcanic arc, suggesting the possible presence of fluids and/or melt at these depths.Mechanisms involving mantle wedge dynamics and grain size effects may contribute to melt generation and magma migration from these reservoirs to the volcanic fronts (15,36).While our model lacks the vertical resolution to speculate on the migration paths of these magmatic systems from the melt reservoir to the volcano, a checkerboard resolution test suggests that the lateral resolution remains good up to approximately 50 to 60 km depth (fig.S3).
Offshore (negative elevation profiles in Fig. 3), the low-rigidity material of the upper-plate contrasts with the mechanically stronger Pacific slab (blue colors in Fig. 2B and 3).The clear distinction between these two separate domains is confirmed by synthetic tests (figs.S4 and S5), and this feature is responsible for the reduction of geodetic residuals offshore compared to the homogeneous model (Fig. 1), with remaining misfit mainly within uncertainties, except for station GTJ4 (table S1).The imaged large-scale variations in shear modulus substantiate the inferences based on forward tests (24,26).Our model's resolution is limited offshore due to the station coverage and trade-offs between the joint parameters (figs.S3 and S5) (19).However, general patterns are confirmed even if the geometry of the Pacific plate's slab interface is not accurately known (figs.S6 and S7).The more compliant forearc structure, with rigidity reductions exceeding 30%, is well correlated with low seismic velocity anomalies in this area (37) and the region of highest seismicity (black dots in Fig. 3).This may indicate mechanically weaker accretionary materials and sediments containing abundant fluids (38), such as mudstones and pelagic clay recovered from ocean drilling expeditions (JFAST) in the frontal prism after the 2011 earthquake (39).
An intriguing observation from Fig. 3, profile C-D, is the segmentation of the seismic activity, as well as the weaker materials, in the forearc.This relatively higher rigidity channel appears to coincide with a region of low seismicity, possibly indicating the control of the upper-plate structure on seismicity patterns.Overall, our results suggest that while geodetic inversions yield consistent estimates of fault slip distributions only on large spatial scales (18), this smoothness does not imply that all information contained in coseismic surface displacements is exploited by standard elastic modeling.Additional insights into Earth's structure and dynamics can be gained by inverting for the medium surrounding the fault.

Implications of geodetic constraints for material properties
Seismic tomography is commonly used to infer lithospheric and mantle structure in subduction zones based on wave speed variations.These anomalies can be further converted into, e.g., temperature or density variations based on a range of assumptions.Shear-wave velocities depend on the square root of the shear modulus divided by density; for instance, the ~30% rigidity reduction beneath the volcanoes that our model recovers (Fig. 4A) corresponds to ≈16% shear-wave reduction for constant density.Variations up to that order are found in the shear-wave anomalies in regional tomographic inversions (Fig. 4B), recognizing that both tomography and rigidity inversions are subject to damping, which leads to underestimations of anomaly strength, and spatially variable coverage and hence resolution.
Our method stands apart from traditional seismic inversions because it offers direct inferences of shear modulus derived from geodetic data, which adds further constraints on the properties of interest.For example, combining our joint inversion with tomographic models in well-constrained regions may allow us to infer density variations more robustly.The density estimate of Fig. 4C is based on combining our inferred 3D rigidity and a smoothed version of shear-wave tomography (33), where smoothing is required to achieve comparable representations since seismic models provide finer-scale features than our geodetic inference.
Negative density anomalies in Fig. 4C appear primarily associated with sediments and weaker rocks offshore, as well as regions where we expect magmatic reservoirs beneath the volcanoes.Positive density anomalies are found in the subducting slab, as expected.The robustness of second-order variations could be further explored in conjunction with gravity constraints (16).More generally, integrating different geophysical information holds the potential to yield more robust estimates of other dynamically relevant parameters for subduction zone dynamics, including temperature and the degree of partial melting.
The complex and highly heterogeneous elastic structure in northeastern Japan we imaged thus complements information from seismic tomography (17,33,40), while providing better constraints on coseismic fault behavior (24,27,(41)(42)(43).This also emphasizes the critical role of material heterogeneity in fault slip inversions and potentially inferring fault stress state, in particular in the presence of a lower-rigidity volcanic arc and upper plate, as well as a stronger slab (24,26,27).Such material property anomalies also appear associated with a heterogeneous viscoelastic, postseismic response where low viscosity beneath the volcanic arc has been invoked to account for subsidence around Quaternary volcanoes (44)(45)(46).
Our joint inversion for 3D elastic variations provides insights into lithospheric dynamics beneath northeastern Japan.We identified low-rigidity material beneath volcanoes, offering a new tool, e.g., for potentially imaging time-dependent magmatic systems.More generally, our approach complements other geophysical inversion techniques.Rigidity imaging has the potential to be integrated with seismic and electromagnetic inversions, enabling a comprehensive and mechanically consistent characterization of thermomechanical subduction zone structure, for example.Our method promises to be especially powerful in settings where densely sampling InSAR deformation time series are available, suggesting the potential for new ways of constraining fault and volcano dynamics.

Cosesmic displacement data
We use the coseismic deformation data compilation of Hashima et al. (24), which consists of 1283 terrestrial geodetic stations managed by the Geospatial Information Authority of Japan.This agency regularly releases daily site location solutions estimated through routine analysis.The effects of the aftershocks were removed following Nishimura et al. (23).In addition to the onshore sensors, we incorporated data from eight acoustic GPS sensors (GPS/A) and six pressure gauges.Among the eight acoustic GPS sensors, six were deployed by the Japan Coast Guard, while the remaining two were installed by Tohoku University (20,21).Tohoku University also deployed four out of the six pressure gauges (22,47), whereas the other two were deployed by the Earthquake Research Institute of the University of Tokyo (48).
One of the GPS/A stations, GJT3, recorded vertical displacements using both the acoustic GPS and a pressure gauge.We opted to use the pressure gauge data for vertical displacement, considering its higher accuracy compared to the GPS/A measurements (24).The offshore stations exhibit larger errors compared to the onshore data (see table S1) (20,24,25).The GPS and APG sensors provide continuous data, while GPS/A observations are from campaign mode measurements (25).Besides primarily recording the Tohoku-oki M9 coseismic displacement, the GPS/A data collected soon but not immediately after the M9, may therefore also contain signals from foreshocks, aftershocks, and postseismic effects (20,21).However, the postseismic contributions should be relatively small (a few tens of centimeters) and thus likely within uncertainties (20,21).Thus, while our results may be affected by undersampled transients like other coseismic study using these data, results are unlikely to be notably biased in this sense (see fig.S3).

Model setup
We constructed a 3D Cartesian model of Japan, which incorporates the Pacific plate boundary and the subducting slab.The fault interface geometry was derived based on interplate seismicity, following the approach outlined in Hashima et al. (24), to allow for comparison.We do not incorporate the Philippine Sea plate and the Nankai Though, as previous studies have demonstrated their minimal impact on the overall coseismic slip distribution during the Tohokuoki earthquake (24).
We converted the spherical coordinates to Cartesian coordinates using an azimuthal equidistant projection with a fixed center point at (140°E, 40°N).For simplicity, we also use a flat surface, disregarding both topography and bathymetry effects.Neglecting sphericity typically induces small changes in surface displacements (49).These errors, comparable to those observed in offshore noisy data (table S1), are not expected to notably alter rigidity variation patterns, although minor amplitude errors may result.Neglecting the effects of topography could potentially lead to bias in slip inversions (50,51), and we expect such effects to be most important offshore where bathymetry changes more strongly than topography on land.Offshore, our model has poor resolution due to the station coverage (fig.S3) and bathymetric gradients might further obscure rigidity structure beyond the trade-offs with fault slip (19).We intend to explore such more example trade-offs in future Bayesian inversions.
On the basis of the fault interface geometry, we generated a 3D mesh using the open-source mesh generation software GMSH (52).The mesh was refined around the Tohoku-oki earthquake epicenter and the central part of eastern Honshu (fig.S1), where we expect most of the information that can be recovered from the surface geodetic data (fig.S3).Our model domain is divided into ≈150,000 tetrahedral elements.The characteristic length of tetrahedra in the mesh ranges from ~10 km near the trench to ~200 km near the lateral and bottom sides of the model domain.For the region of interest, we kept the smallest elements of 10 km fixed up to ~100 km depth to ensure good vertical resolution, gradually increasing toward the bottom of the domain.We have verified that further mesh refinement only yielded negligible changes in the predicted displacement field.
Specifically focusing on the fault, we applied additional refinement near the earthquake location and the areas where we anticipated the highest information recovery concerning slip and subduction zone structure.In these regions, the smallest triangular element on the fault interface had a length of 10 km, gradually increasing to approximately 40 km when moving away from the targeted region.We also confirmed that this choice of reducing the fault mesh resolution outside the expected Tohoku-oki earthquake slip area had no impact on our inverse results.

Forward problem
To solve the elastic forward problem, we follow Puel et al. (53) using a mixed finite element approach that incorporates a fault discontinuity (53), with an implementation based on the open-source FENICS library (54).Our forward model was discretized using a second-order stable triplet of finite-element spaces (19,53), resulting in ≈12.5 million degrees of freedom for the primary variables, namely, stress, displacement, and rotation.For all our tests, we imposed zero displacements on the lateral and bottom boundaries and verified that these boundaries are far enough from the region of interest to affect the model results.
To represent the slip vector s = (s strike , s dip ) on the fault interface, we used the fault local coordinate system instead of expressing the components in the global model coordinates.We restricted the slip to occur only in the strike and dip directions, not allowing opening in the perpendicular direction to the fault plane.To write the slip in its strike and dip directions, we need to substitute the fault term in the elasticity variational form in Puel et al. (53) with where Γ F is the fault internal boundary. and . n × z indicates the cross product between the unit normal vector to the fault plane and the vertical basis z = (0,0,1).τ is a test function, while n ΓF indicates the unit normal to the fault plane, and dS is the integration over the fault boundary.We used the sparse direct solver MUMPS to solve both the forward and inverse problems.

Coseismic slip inversion
We used the open-source library HIPPYLIB (55,56) for all inversions.
Gradient and Hessian information of the least-squares cost functional that penalizes a combination of data misfit and model roughness are computed using adjoint-based optimization methods (57).The fault slip inversion was performed without computing Green's functions following Puel et al. (53).We used a Tikhonov regularization that include a weighted sum of the squared L 2 norm of slip gradient with penalty term γ and the squared L 2 norm of the slip itself with penalty term δ.The first term promotes smoothness of the fault split solution, and the second term ensures strong convexity of the regularization term.To determine suitable penalty weights, we fixed the ratio γ/δ at 10 9 to allow variations of slip of the order of ~30 km.Subsequently, we performed an L-curve analysis (58) to identify the preferred value of γ, which was found to be ∼100.
We scaled the data misfit by the data noise covariance given the observational errors in table S1.We assigned different weights to the horizontal and vertical components of the displacement field to reproduce seafloor deformation, as suggested by several studies (24,25,59).To test the effect of these weights, we performed several coseismic slip inversions assuming a homogeneous medium, computing the L 2 norm of the residual, weighted by the noise uncertainties.By examining the intersection of the different curves, we selected the weighting ratio that minimized this metric.The best ratio was found to be 1: 1  2 : 1 13 : 1 10 for horizontal components of the onshore data, vertical components of the onshore data, seafloor acoustic data, and pressure gauges vertical displacements, respectively (fig.S9).
Using these weighted data, we first performed a coseismic slip inversion for a homogeneous, elastic medium using all available data to determine the fault slip distribution during the Tohoku-oki earthquake, assuming a homogeneous medium with Poisson ratio of 0.25.The slip was discretized using linear Lagrange elements, resulting in approximately 3000 degrees of freedom for each component (strike and dip).We used a conjugate gradient (CG) algorithm preconditioned by the regularization to minimize the gradient (53).The inverse solution converged in ~400 CG iterations, reducing the gradient norm by 12 orders of magnitude.

Joint inversion
The joint inversion method is built upon the approach proposed by Puel et al. (19).For consistency, we used an identical 3D mesh and data weighting as used in the homogeneous coseismic slip problem.Rigidity, represented by the shear modulus parameter of linear elasticity, was determined using a hyperbolic tangent function parameterization to ensure nonnegativity.The resolved patterns remain insensitive to this choice (19).Without stress constraints, our inversion is only sensitive to relative rigidity variations; we constrained the shear modulus to vary between ±50% from a nominal shear modulus background of 60 GPa, while keeping a constant Poisson's ratio of 0.25.
To regularize both the coseismic slip and the shear modulus structure beneath Japan, we used Tikhonov regularization.The penalty weights for the slip components remained the same as in the homogeneous case, while for the shear modulus parameter, we considered a correlation length ( √ γ∕δ ) of ≈30 km, given that the average spacing between land stations is ~20 to 25 km.The value of γ used in our inversion was again determined from L-curve analysis, resulting in a value of 400 (fig.S10).The shear modulus structure was discretized using linear Lagrange elements resulting in ≈25,000 degrees-of-freedom.We used an inexact Newton-CD algorithm (60) to solve the optimization problem, and convergence was achieved within ~20 iterations, reducing the norm of the gradient by six orders of magnitude.
To assess the robustness of our results, we conducted three synthetic recovery and several checkerboard tests.The first one consisted of 3D synthetics where we prescribed the inferred slip from the homogeneous slip inversion and introduced a subduction zone structure with various rigidity values, including a higher rigidity subducting slab (75 GPa), a weaker overriding plate (45 GPa), and lower rigidity spherical anomalies (35 GPa) beneath the volcanic arc within a 30 km radius (fig.S4).We polluted the synthetic data with random Gaussian noise with SD of 1.5 cm.Results show successful recovery of the stronger slab and weaker overriding plate except close to the trench, likely due to trade-offs between slip and material heterogeneity (19).Moreover, our model accurately captured the weaker material beneath the volcanoes in the region of interest up to depths of ≈50 to 60 km.
The effects of the low spatial coverage offshore and trade-offs between the inferred slip and material heterogeneity are more pronounced in a second, antithetical test in which the heterogeneous structure is characterized by a hypothetical weaker subducting slab (35 GPa) and a stronger upper-plate (75 GPa) and volcanic arc (85 GPa) (fig.S5).Although the results show an overall good recovery of the higher rigidity in the overriding plate and the weaker subducting slab in the region of highest resolution (fig.S3), the joint inversion seems to prefer the recovery of a weaker forearc and a stronger slab near the trench.
The third test repeated our reference joint inversion but increased the Pacific plate dip by 10° to evaluate the robustness of the results when the fault geometry is uncertain.The results from this test in figs.S6 and S7 show similar rigidity results as in our reference case (Figs. 2 and 3) with only minor changes in the subduction zone structure and slightly higher slip toward the trench (fig.S2C).
Model resolution was evaluated through a checkerboard test (fig.S3).The test used a sinusoidal model perturbation with a maximum shear modulus variation amplitude of ±35% in a 50 km × 50 km × 50 km grid.On the basis of the recovery, our model provides regional constraints onshore above depths of 40 to 50 km, with good resolution extending to 50 to 60 km beneath the volcanoes.Offshore, resolution is poor due to limited station coverage, with the highest resolution occurring at the intersection point of the profiles (143°E, 38°E) up to depths of 30 to 40 km.

Density calculation
To estimate the density variations from rigidity and shear-wave velocity variations, we take the derivative of the v s formula, v s = √ μ∕ρ , with respect to both shear modulus and density.Expressed as relative variations, the relationship reads: δρ = δμ − 2δv s .

Supplementary Materials
This PDF file includes: Figs.S1 to S10 table S1

Fig. 1 .
Fig. 1.Coseismic geodetic data, displacement residuals, and inferred fault slip distribution for homogeneous and heterogeneous structure.(A) coseismic deformation of the 2011 tohoku-oki M9 earthquake (20-22).Yellow and white arrows indicate horizontal displacement vectors on different scales and uplift is represented by the background color.the red star is the tohoku-oki epicenter, and the purple contour represents the 5-m slip contour from an earlier coseismic slip inversion (24).the dark red lines are major plate boundaries (61), and the red dashed line depicts the 100-km contour line of the subducting Pacific slab (62).(B and C) displacement residuals (on different scale) and fault slip (5 m cyan contours) from our homogeneous medium slip inversion (B), and our joint inversion for fault slip and 3d shear modulus variations (c), as shown in Fig. 2. Mean values in averaging brackets indicate the root mean square (RMS) values (A) and residuals [(B) and (c)] of horizontal (u h ) and vertical (u z ) displacements, restricted to the map area on land.

Fig. 2 .
Fig. 2. Map view of rigidity variations from the 3D joint inversion.the shear modulus, or rigidity, variations (δμ) are shown relative to the background in percent for our preferred joint inversion (Fig. 1c) at depths of 5 km (A) and 75 km (B).Geodetic sites are marked with black circles, volcanoes (63) are denoted by orange triangles.Slip contours in 5-m intervals are represented with black contours.Plate boundaries, slab contours, and epicenter as in Fig. 1.Profiles indicated are shown in Fig. 3.

Fig. 3 . 9 Fig. 4 .
Fig. 3. Vertical profiles illustrating rigidity (shear modulus) variations (δμ) obtained from the 3D joint inversion.the corresponding profile locations are shown on Fig. 2. the orange line represents the topography and bathymetry, while the orange triangles indicate the locations of volcanoes(63).White and gray dots indicate the location of the GPS stations and the APG along the profiles, respectively.the black points correspond to seismicity with a M w > 3.5 from the Japan Meteorological Agency catalogue between 1973 and 2023.the black solid and dashed lines depict the geometry of the Pacific slab and the 100-km depth line, respectively.this depth region represents the level at which the release and upward migration of hydrous minerals' metamorphic reactions within the subducting slab are considered to occur, contributing to the formation of the volcanic arc.